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Abstract 

o 

£T} ■ A systematic and easy-to-use method is developed to calculate directly the doubly heavy 

hadron spectral density in the coordinate space. The correlation function is expressed in terms 
of hypergeometric functions, and the spectral density is obtained through two independent ap- 
proaches: the simple integral representation method and the epsilon-expansion method, respec- 
tively. It is found that the spectral density of doubly heavy hadrons can be analytically expressed 
|- through commonly known simple functions. This method can drastically simplify and improve 

the QCD spectral sum rule calculation of the doubly heavy hadrons. An instructive numerical 
method is also presented for fast evaluation of the spectral density. 
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^ '■ 1 Introduction 

O . 

(f) , In the past years, lots of XYZ states [1-3] have been observed by BaBar and Belle collaborations. 

Many of these states seem not to have a conventional cc or bb structure, which inspired the extensive 
study of the exotic hadron spectroscopy. Then many possible structures like tetraquark, hybrid or 
CS| ' molecular state are proposed and studied [4-12] with QCD spectral sum rules [13-15]. 

These structures can be formally expressed as {QQ'X}, where Q (Q') denotes heavy quark or anti- 
quark and X denotes light quarks and/or gluons. The QCD spectral sum rules require the knowledge 
of the two-point correlation function or its discontinuity, the spectral density, to study many funda- 
mental properties of hadrons. Spectral density is the most labor-intensive part of the QCD spectral 
sum rule calculation. Therefore, it is important and helpful to develop a systematic and easy-to-use 
method to calculate doubly heavy hadron spectral density. 

For a hadron state of light quarks and/or gluons, the correlation function can be easily calculated 
by a Fourier transformation in the coordinate space [16, 17], because the light quark/gluon masses 
are relatively small and the propagators can be approximated by fractional functions. Once the 
heavy quarks are introduced into this kind of states, one need either to use the momentum space 
representation of the heavy quark propagators to keep the masses finite and then evaluate momentum 
integrals [4, 18], or to handle Bessel function related integrals in coordinate space [19]. 

In recent years, the momentum space method [4, 18] has been widely used for mesons [4-6, 8, 20- 
25] and baryons [26] at the leading order in a s . Even so, it is tedious to evaluate two-loop momentum 
integrals with the traditional methods, especially for integrals with tensor structures. What's more, 
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the double integral representations arc time-consuming and might introduce noticeable errors to the 
sum rules. 

In this work, a direct and simple method of analytically calculating the doubly heavy hadron 
spectral density in the coordinate space is developed. In this method, the heavy quark propagator is 
expressed in the form of the modified Bcssel function of the second kind, and the correlation function 
can be expressed by a few generalized hypergeometric functions. The hypergeometric representations 
make the correlation function easily expressed in a compact form and the spectral density easy to 
be worked out analytically. Note that the negative dimensional integration method (NDIM) [27-29] 
or its optimized version, the method of brackets (MB) [30-33], is used to calculate some important 
integrals. 

In Rcfs. [19, 34, 35], the spectral density is calculated in several special cases, where small momen- 
tum expansion q 2 — > 0, large momentum expansion q 2 — > — oo and threshold expansion are used to 
simplify the calculation. Obviously, the spectral density calculated in these special cases is not valid 
for a wide energy region of the QCD spectral sum rules calculation. 

In this work, the spectral density that valid for a wide energy region is calculated in two independent 
approaches: the simple integral representation method and the e-expansion method, respectively. The 
simple (onefold) integral representation [34] of q +\F q type hypergeometric functions is suitable for 
spectral density calculation, where the knowledge of the power function's discontinuity is enough to 
obtain the result. In the recent decade, lots of algorithms or packages [36-42] have been developed to 
perform the e-expansion of hypergeometric functions. Practically, HypExp [38, 39] is used to perform 
e-expansion of q +iF q type hypergeometric functions. It is found that the spectral density of doubly 
heavy hadrons can be analytically expressed in terms of commonly known simple functions and no 
parameter integral is needed at all. 

The well regularized hypergeometric representation of the spectral density can be evaluated nu- 
merically, too. An instructive method is developed for numerical computation of the spectral density. 
In this method, no analytic e-expansion is needed, and one can treat hypergeometric functions as 
common functions. Consequently, the spectral density can be computed directly, which will be helpful 
to the standard community. 

Generally, the coordinate space method can be widely applied to the QCD spectral sum rule 
calculations of doubly heavy hadron states {QQ'X}. This method provides a systematic and easy- 
to-use approach to calculate directly the doubly heavy hadron spectral density in the coordinate 
space, which drastically simplifies the QCD spectral sum rule calculation. Expressing the spectral 
density in terms of simply functions makes the QCD sum rule calculation extremely efficient. There 
is no noticeable errors from multi-dimensional numerical integrals, cither. Then a Monte-Carlo based 
uncertainty analysis [43, 44] is feasible to make realistic uncertainty estimates of the phenomenological 
parameters. As the tetraquark or molecular structure of these states leads to (almost) the same 
predictions within the accuracy [7], such analysis can quantitatively improve the predictive ability of 
QCD sum rules. Moreover, this method can also serve to be an important cross-check of the widely 
used momentum representation method. 

The paper is organized as follows. After a brief introduction of the traditional momentum space 
method in the next section, the method of calculating the doubly heavy hadron spectral density in 
the coordinate space is presented in detail in Sec. 3. The compact hypergeometric representation of 
the correlation function is derived in Sec. 3.2. Then two approaches to extract discontinuities from 
the hypergeometric functions are provided in Sec. 3.3 and Sec. 3.4, respectively. In Sec. 4, a concrete 
calculation is performed to show how to apply the coordinate space method to the practical problems. 
In Sec. 5, an instructive numerical method is also presented for fast evaluation of the spectral density 
from the e-regularized hypergeometric functions. Finally, a summary is given in Sec. 6. 



2 



2 The momentum space method 



The QCD spectral sum rules use the two-point correlation functions like 

iWg) = i [ ^xe^(0\T{j^x)ji(0)}\0) (1) 



to study lots of properties of hadrons. Generally, the correlator in coordinate space can be expressed 
as a product of full propagators and/or their derivatives [19]. And then the correlation function in the 
momentum space can be expressed as a Fourier transformation of the product of several propagators. 
For example, the Lorentz invariant two-point correlation function of {QQqq} state in coordinate space 
is of the simple form 

n(g 2 ) ~ i [ d D xe^S Q (x)S Q (-x)S q (x)S q (-x), (2) 



where all color indices, flavor indices and gamma matrices are ignored, for the sake of briefness. 
The discontinuity of n(g 2 ) is the spectral density, which is needed by the QCD spectral sum rules 
calculation. 

In the past years, the discontinuity of this type of integrals is calculated for mesons [4-6, 8, 20-25] 
and baryons [26] at the leading order in a s by using the techniques in Ref. [4, 18]. To keep the heavy 
quark mass finite, the momentum space expression of the heavy quark propagator [15] 

S<lb i Kb ^"(P + m ) + (P + m )°* V 



4 2 »' v (p 2 -m 2 ) 2 
, i6ab i 2 n 2\ P 2 + mp 

is used, where p = 7 M p M . The light quark part of the correlation function is calculated in the coordinate 
space, and then the whole expression is converted to the momentum space by a D dimensional Fourier 
transformation. Finally, the two-loop momentum integral is calculated and the spectral density is 
expressed as double integrals of Feynman or Schwinger parameters. 
The two-loop momentum integral is of the form 

d 4 k 1 d 4 k 2 jV(k 1: k 2 ,q,m) 

(2tt) 8 [-k 2 + 771 2 ]"i [-jfe| + TO 2 ]"2 [-(q - k~! + k 2 ) 2 }^ ' ( } 

It is tedious to calculate this type of integrals in the traditional manner, and the non-trivial numerator 
Af(ki,k2,q,m) makes the calculation even complicated. Besides, the double integral representations 
of Feynman or Schwinger parameters are time-consuming and might introduce noticeable errors to the 
sum rules. One can also use the Mellin-Barncs method to express such integrals as combinations of 
hypcrgcomctric functions, but the tensor structure makes the final expression diffuse. Therefore, it is 
of great significance to find a convenient way of calculating the spectral density. 

Actually, the spectral densities of doubly heavy hadrons can be calculated in a straightforward way 
in the coordinate space. In the following section, a direct and simple method is developed to calculate 
the doubly heavy spectral density in the coordinate space. 



3 The coordinate space method 

Since the correlation function can be expressed as the Fourier transformation of a product of prop- 
agators in the coordinate space, it is natural to calculate the correlation function in the coordinate 
space. The light quark propagator and the gluon propagator can be expressed as common functions, 
because the masses are small and the small mass expansion can be used. As for the heavy quark, 
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the modified Bcsscl function of the second kind K n (mx) is needed to express the propagator in the 
coordinate space. Note that the so-called fixed-point gauge (x fJ, A )Ji (x) = 0) technique is used to get 
the full propagators. 

3.1 Propagators 

The free-standing part of the light quark propagator in the coordinate space is [16, 17, 44] 

b, 5ab ~ Sab i- S ab g 2 x 2 x 2 



1 ™& + xa» v S ab x 2 , 

9sG ab = + ——(g s qaGq) 



32tt 2; " x 2 ' 192 

6 ab - i 



2 10 x 3 3 



(qq)(giG 2 



m„6 ab m q 5 ab x,_ x m/Vj 4 
" W + —^ {qq) 2^x3* m 

+ I^^V H-x 2 ) - 

-^^^ 2 M-^)(^G 2 ) + ..., (5) 
where x = ix — and G ab v = G™ u T™ b = G™ v \™ b /2. There is also a dangling-gluonic part [17, 44] 



iS%£(x)=Q\T{<f{x)g a (^?m*) 
= - 26 x 3 v»uT2 b {g s q<rGq) 

+ 28 x 9 g + xcr^)T™ b {g s q<jGq) + ■■■ , (6) 

which is important for the condensate contributions from different quarks. 

The heavy quark propagator in the coordinate space can be obtained from the Fourier transfor- 
mation of iSq(p) in (3), 

M$(x) = J ^<(p)e^. (7) 

By using the Fourier transformation relation in D dimensional Minkowski space 

d D p er'P* _ m D - 2v r v - D l 2 K v _ D/2 {r) 
(2ir) D {p 2 - m 2 Y ~ (-2)— 1 {2n) D / 2 T{v) ' 



(8) 



M 



where K v (r) is the modified Bessel function of the second kind, the heavy quark propagator in the 
coordinate space can be expressed as 

iS Q (*) =^ {rr- 2 K 2 (r) + r^K^r)} 

- {{a^r + ra^)r- 1 K 1 (r)+2cr^K (r)} 



576(2tt) 2 to 



{(F-6)r 1 Xi(r)+r 2 ^ 2 (r)} + ••■ , (9) 



where r = mx, r = m\/—x 2 and r 2 = r 2 . 
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In the coordinate space, components of the light quark propagator (5) are of the form (—a; 2 )" or 
x{— x 2 ) n . Note that the (—a; 2 )"" 2 ln(— x 2 ) terms can be rewritten in the form of (— x 2 ) n ~ D / 2 T{D /2 — 
n). As for the heavy quark propagator (9), components are of the form (—x 2 ) n K y (m^f—x 2 ) or 
x{— x 2 ) n K v {m^/ — x 2 ). Naturally, the correlation function (2) can be expressed as D dimensional 
Fourier transformation of two Bcsscl functions. 



3.2 Hypergeometric representation of the correlation function 

After some algebras of the Dirac gamma matrices, the Gcll-Mann matrices and the color indices [45], 
the correlation function reads 

n^g) ~ i /" d D xe**{-gT,^}Y\V :: ^K Vi (m\f^)K Vi (my/^), (10) 
Jm 

where {—g^, x^x"} are possible Lorentz structures and "M" indicates the above integral is calculated 
in the Minkowski space. The angular integral of the Fourier transformation is trivial, 

if d D xe i «-*f(V Z x 1 ) = (2K) D/2 Q 1 - D/2 f°° dxx D / 2 J D/2 _ 1 (Qx)f(x), (11) 
Jm Jo 

where Q 2 = —q 2 , Q > 0, and J v {x) is the Bessel function of the first kind. 

Consequently, the correlation function n(q 2 ) turns into one dimensional integrals of one Bessel J 
and two Bessel K's: 

/>oo 

dx x u ~ x J v (Qx)K a {mx)Kb (mi) . 



This type of integrals can be worked out analytically with the help of negative dimensional integration 
method (NDIM) [27-29] or its optimized version, the method of brackets (MB) [30-33]. In this 
calculation, Bessel functions are needed to be expressed in the form of series. In detail, 

J.(x;) = ^> (-) , (12) 

m— 

where (f> m = jfe^rn ■ Since K v (x) does not have a single summation series representation, one needs 
to express K v (x) definite integral 

1 f°° 

K v (x) = =- dtt^e-iV+V, (13) 

2 Jo 

or as a double summation series 

1 v— v /a;\ ,ll +" 2 
K ^ x ) = 2 {2) (ni-n 2 + v), (14) 

where ^m,n a = ^ni^ an( ^ a formal symbol (a), the bracket, is introduced. This bracket is a short 
form of the divergent integral 

poo 

dxx a - 1 = {a). (15) 



After a little algebra, one gets the key integral of the correlation function calculation in coordinate 
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space 



Here the notation 



dxx u 1 J v (Qx)K a (mx)Kb(mx) 



V T( A__)T(A + _)T(A_ + )T(A ++ ) 
r{l + v)T(u + v) 



X 4^3 



1 I „ u+v l+u+v 
j- t u, 2 ' 2 



-Q 2 \ 

Am 2 J 



(16) 



.4 



AiA 2 



u + v + Aia + A26 



is used to keep the expression short. A similar integral with one Bessel K has been given in the form 
of 2F1 in Ref. [19], which can be used to study hadrons with one heavy quark like D and B mesons 
as well as A c and A& baryons. Since Feynman integrals can be expressed as linear combinations of 
hypcrgcomctric functions, and two-loop self-energy integrals as well as two-loop sunset-type diagrams 
have 4F3 representations [34, 46], the above expression appears as expected. 
Alternatively, the integral above can also be expressed as 



dxx w 1 J v {Qx)K a (mx)K})(mx) 



1 l , r(-a)r(-6)r( 

--2 u ~ 3 m a+b Q~ u - a - b — — — — 



u+v+a-\-b ^ 
2 ' 



_ u—v+a+b 



X4F3 

+ (a -+ -a) + (6 



l+a+b 1 
2 ' x 

1 + a, 1 + 6, l + a + b 



a+p u—v+a+b u+v+a+b 
2 1 2 ' 2 



4m 2 



-b) + (a 



-a, b 



-b). 



(17) 



These two representations (16) and (17) are simply related to each other by analytic continuation of 
the hypergeometric functions. 



q+lFq 



01, . . . , a q +i 
h,...,bg 



r(h) ■ ■ ■ r(b q ) g r(a,o Yltl^i r ^ ~ a *) 



r(oi) • • • r(o g +i) 

x (-z)- a > q+1 F q 



n^r^-a,) 

dh {1 + a » - &k}fc=l,...,g 
{1 + a» — CLk}k=l,...,q+l;k& 



(18) 



where aj—ai ^ Z. In the following calculation, (16) is used to derive the hypcrgcomctric representation 
of the correlation function because this representation can be simply regularized by performing the 
coordinate space integral in D = 4 — 2e dimension. 

Using (11) and (16), the hypcrgcomctric representation of the correlation function can be worked 
out as 



PK =i / d D xe iq - x V-x? K a {m,yj^)K h (m\J^x 2 ) 



M 

_2-D+u)-2 D/2 

r(B__)r(B + _)r(s_ + )r( J B ++ ) 
,B+»»r(f)r(D + ii;) 



x 4^3 



m 

B — , B +A 

D D+w D+w+1 
2 ' 2 ' 2 



(19) 
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where z = s/(Am 2 ) and s = q 2 = —Q 2 . Here the notation 



B 



D + w + X 1 a + X 2 b 



AiA 2 



is used to keep the expression short. Note that D appears in every Gamma function and every 
parameter of the hypergeometric function, which regularizes the divergences very well. It is helpful 
to use the dimensionless variable z as the argument of hypergeometric function. Consequently, all 
dimensional constants like masses and condensates appear in the overall factor. What's more, z 
appears only in the hypergeometric function, which makes it transparent that the correlation function 
has a discontinuity only for z > 1 (or s > Am 2 ). 

In the spectral density calculation, the coordinate space integral is performed in D = 4 — 2e 
dimension, while the propagators of light and heavy quarks in 4 dimension are employed. That is, w, 
a and b in (19) are all integers. The dimension parameter D in the power of constants like 2 D+W ~ 2 , 
-k d I 2 and m D+w can be set to 4 safely, because the e related higher-order infinitesimal terms of the 
constant coefficients will not lead to finite contributions to the spectral density. Furthermore, 4F3 will 
reduce to 3-F2 or even much simpler functions, depending on the specific situation. 

The correlation function contains Lorcntz indices /1 and v, which means Yl^ v {x) contains g^ and 
x^x v structures, as mentioned in (10). Terms proportional to g^ v can be directly calculated through 
(19), and terms proportional to be calculated through the differential form of (19). 

The derivative of hypergeometric function is 



3 F 
dz^ 



{bj} 



n{M 



q+lF q 



{di + 1} 

{bj + 1} 



where II{ai} = 11^ dj and TI{ £»j } 
is 



ITj =1 6j. The double partial derivative of hypergeometric function 



i:) 2 



dq^dq u 



4F3 



{Oi} 

{bj} 



4m" 



q»g v njasK + l)} 
4m 4 U{bj(bj + 1)} 

2m 2 n{bj} 



4F3 



( {a, - 


^2} 


q 2 \ 


\{bj- 


^2} 


Am 2 J 



f i a l - 


HI} 


q 2 \ 


\{b 3 ~ 


Hi} 


Am 2 J 



For the hadron state with J PC = 1 ++ , only the (—g^) part is needed in the spectral density calcu- 
lation. As a result, one gets 



PK 2 =i / d D xe lq - x x^x v ^f^x lW 2 K a {m\/^)K h {my]~x 1 -) 



M 



^2 



D+w -3 g / 2 r(j?__)r(g + _)r(i?_ + )r(i? ++ ) 



x 4F3 



m D+w T(l + f )T(D + w) 
, B+-,B-+, B ++ 



1 



D_ D+w D+w+1 
2 ' 2 ' 2 



(20) 



Note that the factor (—g^ v ) is omitted here for the sake of simplicity. 

By using (19) and (20), one can work out the hypergeometric representations of the correlation 
function IIi(g 2 ). The parameters {a^} and {bj} of q +iF q ({a,i}, {bj}, z) type hypergeometric functions 
are linear functions of space-time dimension D = A — 2e, where e regulates UV and/or IR divergences. 

In fact, it is not too hard to express some integrals as generalized hypergeometric functions, because 
it is trivial, to some extent, to transform these integrals into summations. Nevertheless, if these 
hypergeometric functions can not be easily handled, they are nothing more than short notations and 
the problem remains unsolved. Fortunately, with the help of some new technologies, the q +\F q type 
hypergeometric functions can be effectively handled in many ways. 



7 



To extract the spectral density form the hypcrgeomctric representation of the correlation function, 
two independent approaches are presented below. The simple integral representation method in Sec. 3.3 
is easy-to-use, where the spectral density will be expressed as one-dimensional integrals. In Sec. 3.4, 
the e-expansion method of the hypergeometric functions can even express the spectral density in terms 
of some commonly known functions and no integral is needed. 

3.3 Integral representation of hypergeometric function 

Once the correlation function is expressed as hypergeometric functions like in (19) and (20), one can 
use the simple integral representations [34] of 2-F1 and 3F2 



2F1 



0,1,0,2 



3F2 



h 

01,02,03 
01,62 



r(6i) 



r(o 2 )r(oi - a 2 ) Ji 
r(6 1 )r( 02 



dt 



t ai ~ bl (t, — l) b l- a 2-l 

(t - z) ai 



(21) 



r(a 2 )r(a 3 )r(6i + 6 2 - a 2 - a 3 ) 



< J dt 

x 2F1 



(t - z) Q i 

62 — 02, 62 — 03 

b\ + &2 - a 2 - 03 



1 - t 



(22) 



to calculate the spectral density. It is worth noting that these integral representations are peculiarly 
suitable for the discontinuity calculation. In these representations, z only appears in (t — z)~ ai . 
Generally, the discontinuity of (t — z)~ ai is [19] 



1 , s « (z-t)- ai 8(z-t) 

Disc i - z)- ai = , { , - 1 . 

2ni v ; r(oi)r(l-oi) 

Note that for ai = n — e, where n is a positive integer, the expansion [47, 



(23) 



f L — 4 ml 



m—0 



log m (z) 



(24) 



and its differential forms are used to calculate the discontinuities of the above integrals. The "plus" 
distributions are defined as 



dzf(z) 



/(s)-/(0) , , 
dz g{z). 



(25) 



For the high dimensional condensate parts of the correlation function, -BaiA 2 ma Y be positive. The 
"delta" distribution is more important than the "plus" one, because it will lead to e~ 1 S^(z — t) type 
of contributions to the simple integral. Note that the surface terms of the 5^ (z — t) related integrals 
may contain finite contributions to the sum rules. 

If the hypergeometric functions have already reduced to some simple functions like (z— 1)~ n , where 
n is a positive integer, no integral representation is needed and the discontinuity can be calculated as 
follows 



-LDisc(- 2 r)-" = r J-r -/"-'Hz). 
2m [n — 1)! 



(26) 



3.4 Epsilon expansion of hypergeometric function 

The spectral density can even be analytically worked out from the hypergeometric representation of 
the correlation function. In the recent decade, lots of algorithms and packages [36-42] have been 
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developed to perform the e-cxpansion of hypergeometric functions. Among them, HypExp [38, 39] 
provides a systematic and easy-to-use approach to perform e-expansion of q +\F q type hypergeometric 
functions, whose parameters {cii} and {bj} are of the form n + ae or n + ^ + ae. 

In (19) and (20), the overall factor of Gamma functions has non-zero contribution from order 
0(e~ ni ), while hypergeometric function has non-zero discontinuity from order C(e™ 2 ). For the calcu- 
lations in Sec. 4, it is interesting to see that rt\ = n-i, and then the spectral density calculation can be 
drastically simplified by only calculating e~ ni part of the overall factor and e ni part of the hyperge- 
ometric function. Perhaps this favorable condition does not always exist, but it is worth checking it 
beforehand. 

After e-cxpansion, the hypergeometric functions (or the correlation functions) are expressed in 
terms of commonly known functions and harmonic polylogarithms (HPLs) [49, 50]. The HPLs are of 
the form 



H -\i\l-J^)> (27) 



where • • • denotes indices. Note that the argument i^/~z/(z — 1) of HPLs is not in the well-defined 
interval (0, 1) for z > 1, one cannot naively convert this kind of HPLs to commonly known functions 
for now. 

Empirically, it is preferable to keep the HPLs unchanged before the discontinuities are worked out 
because HPLs are compact functions. Otherwise, one need to calculate the discontinuities of a large 
amount of common functions. The discontinuities are calculated through 

i\ 1 tv tii \ v n(2 + ie) - II(z - ie) 
p(z) = — Discn(z = Inn — . (28 

Z7TZ e->0+ ZlTl 

Particularly, 



= e±iJ , forz>l, (29) 

z ± ls \ z-1 

and now the argument of HPLs becomes — * s [zJ{z — 1) + is or \J z/(z — 1) + ie, respectively. 

Since the HPLs are well defined for arguments in the interval (0, 1) [39, 50, 51], while the real parts 
— y/z/{z — 1) € (—oo,—l) and y/z/Jz — 1) £ (l,oo), these HPLs need to be analytically continued to 
the interval (0, 1) before converting them into commonly known functions. The analytic continuation 
sign S in the small imaginary part ide of the argument is important. In the practical calculation, HPLs 
with argument —Wz/Jz — 1) + ie are analytically continued from (— oo, —1) to (0, 1) by the function 
[50] 

HPLAnalyticContinuation[#, AnalyticContinuationSign -> 1, 
AnalyticContimiationRegion -> minf toml] & 



and HPLs with argument y / z/(z — 1) + ie are analytically continued from (1, oo) to (0,1) by the 
function 

HPLAnalyticContinuation[#, AnalyticContinuationSign -> 1, 
AnalyticContimiationRegion -> onetoinf]& 



After the analytic continuation, the argument of all HPLs becomes y/l — l/z, which is in the interval 
(0, 1) for z > 1. 

Now one can use (28) to calculate the discontinuity of the e-expanded expression. It is known from 
the practical calculations that only a few HPLs (see (43)) are needed to express the spectral density 
in a neat form. Since the argument of all HPLs becomes y/1 — 1/z <G (0, 1) for z > 1, it is safe to 
convert HPLs to commonly known functions like logarithm or polylogarithm. Besides, HPLs can also 
be numerically evaluated with high efficiency like the common functions [50] , and one can directly use 
HPLs in the numerical calculation. 
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The discontinuity near the lower threshold z = 1 is related to functions of the form (1 — z) a . If 
a is a positive integer, one can use (26) to calculate the discontinuity. The balance parameter 

q 9+1 

j=l i=l 

is used to determine the z — > 1 behavior of q +\F q type hypcrgcomctric functions [52. 53]. If a < 0, the 
singular parts of the hypergeometric function are of the form (1 — z) a+k , an additional regularization 
is needed to cancel such kind of divergences. Such contributions might appear when the higher 
dimensional condensates are taken into consideration. From (19) and (20), it is easy to see that 
a = (1 — D)/2 — w or a = (3 — D)/2 — w. Since w is an integer and D will be set to 4 after the 
e-expansion, a cannot be a negative integer. As a result, (1 — z)~ n will not appear in the correlation 
function, and the spectral density will not contain Dirac delta function related terms, either. 



4 Application 

In this section, a concrete calculation is performed to show that the spectral density of doubly heavy 
hadrons like hidden-charm tetraquark state [4] can be easily obtained by using the coordinate space 
method. 

The two-point correlation function (1) in the QCD spectral sum rules can be decomposed into two 
parts as 

rW?) = {^f- - 9^) ni(? 2 ) + ^n ( 9 2 ). (31) 

Tli(q 2 ) and IIo(<? 2 ) have the quantum numbers of the spin 1 and mesons, respectively. For the QCD 
spectral sum rule calculation of JT(3872) meson state with J PC = 1 ++ , only Hi(q 2 ) is of interest. The 
correlation function and the spectral density are related by a dispersion relation 

n 1 ( 9 2 )= [°°ds-^, (32) 
where s< is the lower threshold and the spectral density is the discontinuity of the correlation function 

p(s) = Disc III (s). (33) 
The current of the 1 ++ tetraquark [qc] [qc] is given by [4] 

JV = .fa b ;de{l T aCT A C b ){q d T B Cc T e ), (34) 

where f a b-de is the color structure factor, and the Dirac gamma matrix T is of the form /, 7^, g^ v , 
7 5 7 M or i"f 5 . Note that T's satisfy the transformation relation 

7°r 7 ° = r t . (35) 

The vacuum expectation value of the currents is 

(0|T{j„(x)jt(0)}|0) =<0|T{j„(x)it( )}|0) 1 

+ (r A o r B ) + (r c o r D ) + (r A o r B , r c o r D ). (36) 
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The labeled gamma matrices T A , T B , T c and T D are introduced to make the sources of the gamma 
matrices explicit. In the final expression, these gamma matrices will be replaced by 275, 7^, 175 and 
7„, respectively, for the specific spectral density calculation in Ref. [4]. 

After contracting the quark fields in the correlation function with Wick theorem, one part of the 
vacuum expectation values reads 

(0|r{j M (x)jt(0)}|0)! =(0| Tr{T A [ l S^(x)}T c (C[ l Sl,(x)fC- 1 )} 

x Tr{r B (C[ i 5 e c , e (- 2 ;)] T C- 1 )r^[z^ d (-.T)]}|0), (37) 

where iS^ a , (x) is the full quark propagator in coordinate space. It is easy to obtain the other three 
parts by permutations of gamma matrices. 

In the propagators (5) and (9), the gamma matrix related structures are 1, x, a^ v ' . a^x + xa^ v . 
Their charge conjugation transformations are 

C{1, x, a^,a^x + xa^YC' 1 = {1, —x, -a^,a^x + xa^}. (38) 

With the help of this identity, it is easy to get the charge conjugation of the quark propagators. 
After some algebras of the gamma matrices, the Gell-Mann matrices and the color indices [45], the 
correlation function can be obtained easily. For example, the perturbative part of the correlator in 
coordinate space is 

nf rt ^» = ^ (- 5 ^(-.t 2 )- 4 - 2xW(-x 2 )- 5 ) A^ 2 (m cV /^) 2 . (39) 

By using (19) and (20) with w = —8 and a = b = —2, one gets the II^ olt (2) in the form of hypergeo- 
metric functions. Explicitly, 



3m 



p ( _ e _4)r(- £ -2) 2 r(- e ) /- e -4,-e-2,-e 



Hl (z) -64^\ r(2-£)r(-2e-4) ^ { -e-§,2 



e 



r(-e-4)r(-e-2) 2 r(-e) ( -e - 4, -e - 2, -e 

if-. 



r(3 - e)r(-2e - 4) -e-|,3-e 



(40) 



With the help of (22) and (23), one gets the simple integral representation of the perturbative part of 
the spectral density 



64 ;j[**in«-l) U/ ^l(u 7 



1 - t 



l-t . (41) 



45045 

Here a\ = —e is chosen, which makes the (z — t)~ ai trivial when e expansion is performed. It is 
worth noting that one has the freedom to chose a parameter from {a^} as a\, for {aj} in p F q are 
totally symmetric. With other choices of <xi, different but equivalent integral representations will be 
obtained. Analogously, the other parts of the spectral density can be expressed as simple integrals, 
too. 

Furthermore, the spectral density can also be worked out analytically from (40) by using the e- 
expansion method in Sec. 3.4. Note that the e-expansion of the Gamma functions in (40) begins from 
g_3e -3 , while the non- vanishing discontinuity of the hypergeometric function 3F2 begins from fz(z)t : . 
Technically, it is sufficient to take 5-3/3(2) as the perturbative part of the correlation function, which 
makes the spectral density calculation drastically simplified. Explicitly, the perturbative part of the 
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spectral density is 

P pcrt W 

where 



2 10 7T 6 



32z 4 



48z 2 



1456z 3 776z 2 2624z 7 7 

^ 18 + 12z 



45 
128z 



15 
7 



24z 2 



45 

17(«) + 12T(z) }, 



V(z) 



(42) 

(43) 

(44) 
(45) 



V(z) = y/l- 1/z, U(z) = H + (V(z)), T{z) = H- t+ (V{z)). 
The HPLs [49, 50] H + (z) and H_ t+ (z) can be converted to commonly known functions 

H+{V{z)) =log(l + y(z)) - log(l - V(z)), 

H. t+ (V(z)) =Li 2 (i^W) - Li 2 (±±^>) + | log(4z)t/(z), 

where Li2(z) is the Euler dilogarithm. 

In Rcf . [4] , spectral density is calculated by the momentum integral method. The perturbativc part 
of the spectral density is expressed by double integral of modified Fcynman or Schwinger parameters 

P PeTt ( s ) = oTir~fi / > da I > W-LrO. " « " ^)(! + a + P)K a + P) m l ~ ^l 4 ' ( 46 ) 



a 3 /? 



where a < = (1 — w)/2, a> = (1 + w)/2, /3< = am^/(sa — ml) and /3> = 1 — a, with u = \/l — 4m^/s. 
This double integral representation is obtained by calculating the discontinuity of two-loop sunset type 
momentum integrals like (4). 

V(z) and U{z) vary drastically near the threshold z = 1, which means the numerical computation 
of the double integral representation (46) might contain large relative errors near the threshold. In 
the QCD sum rules approach, high energy part of the spectral density p(z = s/(4m^)) is damped by 
the factor e~ ST , and the large relative errors of the low energy part of the spectral density might cause 
some noticeable errors to the sum rules. The double integrals are also time-consuming. Actually, 
three-fold numerical integrals are involved in the sum rule calculation, if the s integral is taken into 
account. Numerically, (46) is about hundreds of times slower than (42). 

To make realistic uncertainty estimates of the results, a Monte-Carlo based uncertainty analysis 
[43, 44] is often used in the QCD sum rule calculation. In this procedure, the entire phase space of QCD 
input parameters like quark masses and condensates, is explored simultaneously, and is mapped into 
uncertainties in the phcnomcnological parameters. Since a Monte-Carlo based uncertainty analysis 
need more than 100 times of evaluation to get stable uncertainties, it is not wise to repeatedly calculate 
the multi-dimensional integrals. Technically, one can write the double integrals in a dimensionless 
form by taking the dimensional parameters as the coefficients of these integrals. Subsequently, the 
double integrals can be numerically evaluated only once for a sequence of uniformly-spaced points 
z% = Si /(4m;?), and then the Newton-Cotes formula is used to perform the z integration. 

Comparing (41) and (42) with (46), one can find that the coordinate space method is more suitable 
for spectral density and sum rule calculation, at least for the case of two heavy quarks with equal 
masses. 

In the same way, the explicit analytic expressions of the other parts of the spectral density [4] can 
be obtained as 

m q ml J |"48z 3 6172 2 2087z 127 
2 8 tt 6 \ [~ 5~ + "T5 30 2A 

7 



28z' - AOz - 3 + — 



30 
5 



7 



V(z) 



32z 2 



U(z) - 18T(z) 



m q m*(qq) 

2 5 7T 4 



2z l 



37z 

IT 



47 
12 



1 

8z 



9-H 

z 



1 



16z 2 



U(z) 



(47) 
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p™(z) 



p^\z) 



2z 48z 2 v ' 



m 4 c (g 2 G 2 ) 

2 9 3tt 6 
ml(qgaGq) ( 
2 6 tt 4 \ 
m 2 c (qq) 2 

12tt 2 



'4£2 

58z 
~~9~ 



bz 

T ~ : 

7 7 

36 ~ 2Az 



V(z) + 
V(z) 



2z-l 

~Az 3 

y + 2 



1 

Az 



U(z\ 
7 

48z 2 



(49) 
(50) 
(51) 



Accordingly, the spectral densities of the doubly heavy meson state are analytically worked out 
through the coordinate space method. Since no multi-dimensional integral like (46) is involved in 
the spectral density and the analytic functions can be numerically calculated with high efficiency, the 
errors of sum rules from spectral density functions are minimized. Moreover, a Monte-Carlo based 
uncertainty analysis can be performed to make realistic uncertainty estimates of the phenomenological 
parameters, which will improve the predictive ability of QCD sum rules. 



5 Numerical evaluation of spectral density 

As the hypergeometric functions might not be well known for the standard community and the analytic 
reduction might be a little involved, a numerical method is presented in this section for fast evaluation 
of the spectral density from the e-rcgularized hypergeometric functions. Consequently, no analytic 
e-expansion is needed, and the regularization scheme is no longer restricted by the capabilities of 
HypExp [39]. 

Since the spectral density itself is finite [54] and independent of the regularization parameter, an 
unorthodox but reasonable regularization scheme can be used to obtain the spectral density. Formally, 
the well regularized spectral density can be expressed in the form of 

m 

p(z,e) = ]Te-"-O + p(z) + 0(e), (52) 

n=l 

where p(z,e) is the approximation of order 0(e) of the true spectral density p(z). If e is numerically 
small enough, p(z,e) can be taken as p{z). It is worth noting that multiple precision computation 
is needed to get rid of roundoff errors from the Y^=i e_ ™ ' part. Roughly, —(m + 1) lg(e) digits of 
precision is enough to obtain the correct results. 

Practically, (19) and (20) with the regularization parameter D = 4 — 2e are not good for numerical 
computation, because q +iF q ({a,i}, {bj}, z) is evaluated through its analytic continuation (18) for z > 1. 
If the parameters a,j — a% G Z, the arguments of gamma functions as well as the parameters of 
q+iFq^a'j}, {b'j}, 1/z) may contain nonpositive integers, and some spurious divergences may arise. 
Then, additional auxiliary parameters are needed to regularize the expression, which may lead to a 
complex calculation. 

Using (18), (19) and (20), it is straightforward to get the hypergeometric representation of the 
correlation function with argument 1/z. In this section, this representation will be used for numerical 
computation. For z > 1, the hypergeometric functions with argument 1/z have no branch cuts, and 
the discontinuities come only from the (— z) a part. With the help of (23), one gets the discontinuities 
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of the integrals of {—g^ v ) and x^x v structures. Explicitly, 

DK =-!- Disci / d D xe lq - x \T^ K a (m^f^)K b (m^f^) 



M 



X 4^3 



l-\-a-\-b _|_ a-\-b ^ j_ w-\-a-\-b D-\-w~\-a J rb 



l 

+ (a^ -a) + {b^ -b) + (a-> -a, -6), (53) 



2 ' ' 2 ' ' 2 ' 2 

1 + o, 1 + 6, 1 + a + 6 



and 



DK 



— Disci / d D xe* q - x x^x''y/^ U ' ^^m^^fmV^) 



M 



I, q n/o n £>+tu+a+b TV-a)r(-6) 



X 4^3 



D-\-w-\-a-\-b i£+a+b^ 

l+a-j-b ^ _L Q +fr m+a+fr .D+tu+a+fr 



+ (a -> -a) + (6 ->• -6) + (a ->• -a, 6 ->• -6). (54) 



2 ' 1 2 ' 2 ' 2 

1 + a, 1 + 6, l + a + 6 



Note that the factor (—g^ v ) is omitted here for the sake of simplicity. These two functions DKq(-D, w, a, 6, m, z) 
and DK2(-D, w, a, 6, m, z) can be well regularized by an unorthodox rcgularization scheme: D = 4, 
a = ao + ae and 6 = 6o + /3e, with a + 1 0, /3 ^= and a ± /3 7^ 0. It is preferable to use irrational a and 
/3 for practical calculation. 

Accordingly, it is easy to write out the spectral density from DKo(-D, w, a, 6, to, z) and DK2(-D, w, a, 6, to, z). 
It is better to show the utility of this algorithm by a concrete example. From (39), one can write out 
the spectral density 

3m 4 

p pcit (z,e) = —±{DK (D,w,a,b,m c ,z) - 2DK 2 (L>, w, a, 6, m c , z)} (55) 

with D = 4, w = —8, a = —2 + e and 6 = —2 + 3e. Due to the T(—n + ae) terms, a high precision 
computation is needed to get rid of roundoff errors. Fortunately, the hypergeometric function can be 
evaluated to arbitrary numerical precision by computation systems like the Python library mpmath [55] 
or Mathematica. Specifically, some numerical results are listed in Table 1 . Note that the charm quark 
mass to c is set to 1, and the working precision is set to 32 digits for e = 10~ 8 . The exact value of 
p pert (z) can be evaluate by using (42) or (46). 



Table 1: Numerical evaluation of spectral density p pcrt (z) with m c = 1. 





z = 1.1 
(xl(T 13 ) 


z = 2 
(xl(T 7 ) 


z = 5 
(xlO- 4 ) 


z = 10 
(xlO" 2 ) 


z = 20 
(xlO- 1 ) 


e = 1(T 4 
e = 1(T 6 
e = 10- 8 
Exact 


3.847514018 
3.847553295 
3.847553688 
3.847553692 


3.484073813 
3.484354676 
3.484357485 
3.484357513 


3.790001592 
3.790789489 
3.790797370 
3.790797449 


1.489432452 
1.489904972 
1.489909698 
1.489909746 


3.687774818 
3.689380738 
3.689396801 
3.689396963 



From Table 1 , it is easy to see that one can use a finite small parameter to regularize the spectral 
density and obtain the numerical result to a desired precision. Several different input values of e 
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parameter are used to show the numerical stability of the algorithm. As e decreases to zero, the 
numerical result steadily approaches the exact value, yet more digits of working precision is needed. 
Therefore, it is preferable to use a fairly small e to save the computation time. 

The numerical method can also be applied to the calculations of more complicated hypcrgcometric 
functions. 

6 Conclusion 

In this paper, a systematic and easy-to-use method is developed to calculate directly the doubly heavy 
hadron spectral density in the coordinate space. The correlation function is expressed in terms of 
generalized hypergeometric functions, and then the spectral density is obtained through two indepen- 
dent approaches: the simple integral representation method and the e-expansion method, respectively. 
It is found that the spectral density of doubly heavy hadrons can be analytically expressed through 
commonly known functions. After that, a concrete calculation is performed to show that the spec- 
tral density of doubly heavy hadrons can be easily obtained by using the coordinate space method. 
This method can drastically simplify the QCD spectral sum rule calculation of the {ccX} and {bbX} 
systems. 

An instructive numerical method is developed to evaluate the e-rcgularized spectral density directly. 
This method can also be applied to other calculations where expressions are in the form of regularized 
(hypergeometric) functions. In particular, the numerical computation can be used to handle problems 
that are beyond the reach of the analytic e-expansion method. 

Since the spectral density is expressed in terms of simple functions, there is no noticeable errors 
from multi-dimensional numerical integrals. A Monte-Carlo based uncertainty analysis [43, 44] can be 
used to make realistic uncertainty estimates of the phenomenological parameters. Such analysis can 
quantitatively improve the predictive ability of QCD sum rules. Moreover, this method can also serve 
to be an important cross-check of the widely used momentum representation method. 

Furthermore, for a hadron state with two different massive quarks, the correlation function can be 
expressed in the form of Appell F4 functions of two variables, as is shown in A. This type of spectral 
density can be used to study the f? c -like {cbX} structures [21, 25, 56]. It is of great worth to work out 
the spectral density from the Appell F4 functions. The calculations will be more complicated because 
the Appell F4 is intricate and not widely studied like the generalized hypergeometric functions p F q . 
The method can also be extended to {QQQX} and {QQQ'X} systems. These issues will be discussed 
in detail in a future publication. 

Acknowledgments 

Z.W.H. thanks Ze-kun Guo for helpful discussions. This work is supported by the National Natural 
Science Foundation of China under Grant No. 10775105, BEPC National Laboratory Project R&D 
and BES Collaboration Research Foundation, and the project of Wuhan University of China under 
the Grant No. 201103013 and 9yw201115. 

A Two massive quarks with different masses 

In the case of two heavy quarks with different masses, one can also use the method of brackets [30-33] 
to calculate the integral of Bessel functions. Explicitly, this type of integral can be expressed in the 
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form of Appell F4. 



dxx u 1 J v (Qx)K a (mix)Kb(m2x) 
=^-> nh ---^m a -^- (b)r(Af - )r(A -- ) 



+ (b^-b). 



r(i 



f 4 





-Q 2 


4) 


VI 




m{) 



(56) 



Appell F4 is defined as 



Fa 



a; (3 
7i,72 



E 



(/?) 

x y . 



m!n!(7i) m (7 2 ) n 



(57) 



The above integral can also be expressed as Appell F4 with arguments (— Q 2 /m|, mf/mf) or 
(rn\/{— Q 2 ), m\j (— Q 2 )). These three representations are related to each other by analytic continua- 
tion. If mi = 7712 = m, F4 turns into 4-F3. 

With the help of XSummer [57] and/or HYPERDIRE [42], F4 may be reduced to commonly known 
functions and HPLs. If the analytic e-expansion procedure fails, one can use the method in Sec. 5 to 
evaluate the spectral density numerically. 
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